/* Bessel J_1(x) function in double precision */

#include <math.h>

double dbesj1(double x)
{
    int k;
    double w, t, y, v, theta;
    static double a[8] = {
        -1.4810349e-13, 3.363594618e-11, 
        -5.65140051697e-9, 6.7816840144764e-7, 
        -5.425347222188379e-5, 0.00260416666666662438, 
        -0.06249999999999999799, 0.49999999999999999998
    };
    static double b[65] = {
        2.43721316e-12, -9.400554763e-11, 
        3.0605338998e-9, -8.287270492518e-8, 
        1.83020515991344e-6, -3.219783841164382e-5, 
        4.3795830161515318e-4, -0.00442952351530868999, 
        0.03157908273375945955, -0.14682160488052520107, 
        0.39309619054093640008, -0.4795280821510107028, 
        0.1414899934402712514, 
        1.82119257e-12, -6.862117678e-11, 
        2.1732790836e-9, -5.69359291782e-8, 
        1.20771046483277e-6, -2.020151799736374e-5, 
        2.5745933218048448e-4, -0.00238514907946126334, 
        0.01499220060892984289, -0.05707238494868888345, 
        0.10375225210588234727, -0.02721551202427354117, 
        -0.06420643306727498985, 
        1.352611196e-12, -4.9706947875e-11, 
        1.527944986332e-9, -3.8602878823401e-8, 
        7.82618036237845e-7, -1.23499947484511e-5, 
        1.45508295194426686e-4, -0.001203649737425854162, 
        0.006299092495799005109, -0.016449840761170764763, 
        0.002106328565019748701, 0.05852741000686073465, 
        -0.031896615709705053191, 
        9.97982124e-13, -3.5702556073e-11, 
        1.062332772617e-9, -2.5779624221725e-8, 
        4.96382962683556e-7, -7.310776625173004e-6, 
        7.8028107569541842e-5, -5.50624088538081113e-4, 
        0.002081442840335570371, -7.71292652260286633e-4, 
        -0.019541271866742634199, 0.033361194224480445382, 
        0.017516628654559387164, 
        7.31050661e-13, -2.5404499912e-11, 
        7.29360079088e-10, -1.6915375004937e-8, 
        3.06748319652546e-7, -4.151324014331739e-6, 
        3.8793392054271497e-5, -2.11180556924525773e-4, 
        2.74577195102593786e-4, 0.003378676555289966782, 
        -0.013842821799754920148, -0.002041834048574905921, 
        0.032167266073736023299
    };
    static double c[70] = {
        -1.185964494e-11, 3.9110295657e-10, 
        1.80385519493e-9, -5.575391345723e-8, 
        -1.8635897017174e-7, 5.42738239401869e-6, 
        1.181490114244279e-5, -3.300031939852107e-4, 
        -3.7717832892725053e-4, 0.01070685852970608288, 
        0.00356629346707622489, -0.13524776185998074716, 
        0.00980725611657523952, 0.27312196367405374425, 
        -3.029591097e-11, 9.259293559e-11, 
        4.96321971223e-9, -1.518137078639e-8, 
        -5.7045127595547e-7, 1.71237271302072e-6, 
        4.271400348035384e-5, -1.2152454198713258e-4, 
        -0.00184155714921474963, 0.00462994691003219055, 
        0.03671737063840232452, -0.06863857568599167175, 
        -0.21090395092505707655, 0.16126443075752985095, 
        -2.19760208e-11, -2.7659100729e-10, 
        3.74295124827e-9, 3.684765777023e-8, 
        -4.5072801091574e-7, -3.27941630669276e-6, 
        3.5713715545163e-5, 1.7664005411843533e-4, 
        -0.00165119297594774104, -0.00485925381792986774, 
        0.03593306985381680131, 0.04997877588191962563, 
        -0.22913866929783936544, -0.07885001422733148814, 
        5.16292316e-12, -3.9445956763e-10, 
        -6.6220021263e-10, 5.511286218639e-8, 
        5.01257940078e-8, -5.22111059203425e-6, 
        -1.34311394455105e-6, 3.0612891890766805e-4, 
        -7.103391195326182e-5, -0.00949316714311443491, 
        0.00455036998246516948, 0.11540391585989614784, 
        -0.04779493761902840455, -0.2283786206653234746, 
        2.697817493e-11, -1.6633326949e-10, 
        -4.3313486035e-9, 2.508404686362e-8, 
        4.8528284780984e-7, -2.58267851112118e-6, 
        -3.521049080466759e-5, 1.6566324273339952e-4, 
        0.00146474737522491617, -0.00565140892697147306, 
        -0.028338820556793004, 0.07580744376982855057, 
        0.16012275906960187978, -0.16548380461475971845
    };
    static double d[52] = {
        -1.272346002224188092e-14, 3.370464692346669075e-13, 
        -1.144940314335484869e-11, 6.863141561083429745e-10, 
        -9.491933932960924159e-8, 5.301676561445687562e-5, 
        0.162867503967639974, -3.652982212914147794e-13, 
        1.151126750560028914e-11, -5.165585095674343486e-10, 
        4.657991250060549892e-8, -1.186794704692706504e-5, 
        0.01562499999999994026, 
        -8.713069680903981555e-15, 3.140780373478474935e-13, 
        -1.139089186076256597e-11, 6.862299023338785566e-10, 
        -9.491926788274594674e-8, 5.301676558106268323e-5, 
        0.162867503967646622, -2.792555727162752006e-13, 
        1.108650207651756807e-11, -5.156745588549830981e-10, 
        4.657894859077370979e-8, -1.186794650130550256e-5, 
        0.01562499999987299901, 
        -6.304859171204770696e-15, 2.857249044208791652e-13, 
        -1.124956921556753188e-11, 6.858482894906716661e-10, 
        -9.49186795351689846e-8, 5.301676509057781574e-5, 
        0.1628675039678191167, -2.185193490132496053e-13, 
        1.048820673697426074e-11, -5.132819367467680132e-10, 
        4.65740943737299422e-8, -1.186794150862988921e-5, 
        0.01562499999779270706, 
        -4.74041720979200985e-15, 2.578715253644144182e-13, 
        -1.104148898414138857e-11, 6.850134201626289183e-10, 
        -9.49167823417491964e-8, 5.301676277588728159e-5, 
        0.1628675039690033136, -1.75512205749384229e-13, 
        9.848723331445182397e-12, -5.094535425482245697e-10, 
        4.656255982268609304e-8, -1.186792402114394891e-5, 
        0.01562499998712198636
    };

    w = fabs(x);
    if (w < 1) {
        t = w * w;
        y = (((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7]) * w;
    } else if (w < 8.5) {
        t = w * w * 0.0625;
        k = (int) t;
        t -= k + 0.5;
        k *= 13;
        y = ((((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12]) * w;
    } else if (w < 12.5) {
        k = (int) w;
        t = w - (k + 0.5);
        k = 14 * (k - 8);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / w;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= sin(w + theta);
    }
    return x < 0 ? -y : y;
}

